Dynamical phase transition for current statistics in a simple driven diffusive system 
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We consider fluctuations of the time-averaged current in the one-dimensional weakly-asymmetric 
exclusion process on a ring. The optimal density profile which sustains a given fluctuation exhibits 
an instability for low enough currents, where it becomes time-dependent. This instability corre- 
sponds to a dynamical phase transition in the system fluctuation behavior: while typical current 
fluctuations result from the sum of weakly-correlated local events and are still associated with the 
flat, steady-state density profile, for currents below a critical threshold the system self-organizes 
into a macroscopic jammed state in the form of a coherent traveling wave, that hinders transport of 
particles and thus facilitates a time-averaged current fiuctuation well below the average current. We 
analyze in detail this phenomenon using advanced Monte Carlo simulations, and work out macro- 
scopic fluctuation theory predictions, finding very good agreement in all cases. In particular, we 
study not only the current large deviation function, but also the critical current threshold, the asso- 
ciated optimal density profiles and the traveling wave velocity, analyzing in depth finite-size effects 
and hence providing a detailed characterization of the dynamical transition. 
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I. INTRODUCTION 

Recent years are witnessing a quiet revolution in 
nonequilibrium statistical physics. At the core of this 
revolution is the realization of the essential role played 
by macroscopic fluctuations to understand the nonequi- 
librium behavior of a system of interest [IHIO] . This ac- 
tivity has led to a number of groundbreaking results valid 
arbitrarily far from equilibrium (and therefore not re- 
stricted to the confining world of linear response) , which 
are offering a glimpse of the long-sought general theory 
of nonequilibrium phenomena. A main example is the 
Gallavotti-Cohen fluctuation theorem [IHlj, which ex- 
presses the subtle but enduring consequences of micro- 
scopic time reversibility at the macroscopic level. The 
list continues however, with further breakthroughs rang- 
ing from the Jarzynski equality [5 or the Crooks fluctu- 
ation theorem [6] to the Hatano-Sasa relation 7J or the 
recent extension of Clausius inequality to nonequilibrium 
steady states [5], to mention just a few [S]. In addition, a 
general theoretical framework, the macroscopic fluctua- 
tion theory of Bertini and coworkers [10.^, has been devel- 
oped to understand the fluctuating behavior of diffusive 
systems far from equilibrium (with recent generalizations 
to driven dissipative media [TTIIT^ ). 

A general observation underlying many of these re- 
sults is that macroscopic fluctuations are often associ- 
ated with a nontrivial and well-defined path in phase 
space, a path that the system traverses in order to fa- 
cilitate such fluctuation. The properties of these optimal 
paths are revealing a whole new phenomenology at the 
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fluctuating level with important implications out of equi- 
librium |13H15j . For instance, the optimal path leading 
to a macroscopic fluctuation in a nonequilibrium steady 
state has been recently shown to be the time-reversal of 
the relaxation path from this fluctuation according to 
some adjoint hydrodynamic laws (which are not neces- 
sarily equal to the forward-in-time hydrodynamics) |10) . 
This general result valid arbitrarily far from equilibrium 
reduces to the well-known Onsager's regression hypoth- 
esis when small deviations from equilibrium are consid- 
ered. Moreover, the study of the symmetry properties 
of the optimal paths for current fluctuations has led to 
another remarkable insight, the isometric fluctuation re- 
lation [T4j, which in turn implies a set of hierarchies of 
equations for the current cumulants and the nonlinear re- 
sponse coefflcients, going far beyond Onsager reciprocity 
relations and Green-Kubo formulas. 

Another recent and striking discovery concerns the ex- 
istence of coherent structures associated to large, rare 
fluctuations [161 I17j . which in turn imply that these 
events are far more probable than previously anticipated. 
Such coherent, self-organized patterns emerge via a dy- 
namical phase transition at the fluctuating level, which is 
accompanied by spontaneous symmetry breaking |161I17) . 
The aim of this paper is to investigate in detail this phe- 
nomenon in a simple diffusive system in one dimension, 
namely the weakly-asymmetric simple exclusion process 
(WASEP) [IH], where we study fluctuations of the time- 
averaged current. 

The model is defined on a one-dimensional (Id) lattice 
of size N with periodic boundary conditions (pbc), where 
M < N particles live, see Fig. [ija, so the total density 
is po = M/N. Each lattice site i G [1, A^] niay contain at 
most one particle, so the state of the system is defined by 
a set of occupation numbers, n = {ui = 0,1, i £ [1, N]}, 
and M = X)i=i '^j- Dynamics is stochastics and proceeds 
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in nonequilibrium physics a role akin to the free energy of 
equihbrium systems. Therefore, even though we do not 
know how to connect in general microscopic dynamics to 
macroscopic properties in nonequilibrium systems (in a 
way equivalent to the equilibrium ensemble formalism), 
we can still measure LDFs of macroscopic observables 
out of equilibrium, which provide an alternative path to a 
detailed macroscopic description of nonequilibrium phe- 



FIG. 1: (Color online) (a) Sketch of the WASEP. Particles 
in a periodic Id lattice jump stochastically to a right (left) 
empty nearest neighbor at a rate r+ (r-), so particles feel 
an external driving field E = ^ln(^). (b) Convergence of 
the time-averaged current to its ensemble value (g) for many 
different realizations (top line cloud), and sketch of the prob- 
ability concentration as time increases, associated with the 
large deviation principle, eq. ([T]). 



via sequential particle jumps to nearest neighbor sites, 
provided these are empty, at a rate r± = ^ exp{zLE / N) 
for jumps along the ±i;-direction [19]. Here E plays the 
role of a weak external field which drives the system to a 
nonequilibrium steady state characterized by a homoge- 
neous density profile {p{x)) = po and a nonzero net av- 
erage current (q) = po(l — Po)E. We employ continuous- 
time Markov dynamics, so the time to exit a configura- 
tion n is a random variable drawn from a Poisson distri- 
bution with exit rate R{n) = M+r+ -|- M_r_, where M± 
is the number of particles with empty nearest neighbor 
in the ±i-direction. 

We are interested in the statistics of the total parti- 
cle current q flowing through the system, averaged over 
a long diffusive time r and across space. In partic- 
ular, we define the empirical time-averaged current as 
q ~ t^^{Q^ — Q^), where Qf is the total number of 
particle jumps in the ii-direction in a given microscopic 
time interval t ~ tN"^ . For r — ^ oo this time-averaged es- 
timate converges toward the ensemble average {q) . How- 
ever, for long but finite times r we observe fluctuations 
q^ (q), and their probability P,- (g) obeys a large devia- 
tion principle in this limit ^2D] 



Priq) 



(1) 



where G{q) < is the current large deviation function 
(LDF), such that G{{q)) = 0. This means that the prob- 
ability of observing a fixed current fiuctuation q ^ (q) 
decays exponentially as both r and N increase, at a rate 
given by G{q), see Fig. [l]b. In other words, G{q) mea- 
sures the rate at which Priq) concentrates around (q). 
The above large-deviation principle describes the scaling 
oiPriq) for both typical and rare current fluctuations. In 
particular, a suitable expansion of G{q) for small fluctua- 
tions yields the usual gaussian form for Pt('?) associated 
with the central limit theorem. The current LDF plays an 
important role in nonequilibrium statistical physics as it 
contains essential information on the transport properties 
of the system at hand. Moreover, in general LDFs play 



II. MACROSCOPIC FLUCTUATION THEORY 
AND DYNAMICAL PHASE TRANSITION 

Computing LDFs from scratch, starting from micro- 
scopic dynamics, is a humongous task which has been 
achieved only in a handful of oversimplifled models (most 
of them stochastic lattice gases and related models) 
[TUl [TB] . However, in a recent series of works, Bertini 
and collaborators [10] have developed a phenomenolog- 
ical theory, the macroscopic fluctuation theory (MFT), 
which describes in detail dynamic fluctuations in driven 
diffusive systems starting from the hydrodynamic evolu- 
tion equation for the local density p{x, t) for the system 
of interest and the sole knowledge of two transport coeffi- 
cients, the diffusivity D{p) and the mobility a{p), which 
can be measured experimentally. From this knowledge, 
MFT offers explicit predictions for the current LDF (see 
the Appendix) 

Giq) = -i min r dt r ,J^ + md.P-^ip)E? 

(2) 

where a long-time limit is implicit and the minimum is 
taken over all histories of the density and current fields, 
p{x, t) and j(x, t) respectively, coupled via the continuity 
equation dtp + dxj = at every point of space and time, 
and subject to the constraint q = dt dxj{x,t) 

for the space&time-averaged current and the appropri- 
ate boundary conditions (periodic in this case), see Ap- 
"for details. Note that, for the WASEP, D{p) = \ 
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and a{p) = p{\ — p) [TH]. The optimal density and cur- 
rent fields solution of the above variational problem, de- 
noted here as pq{x,t) and jq{x,t), can be interpreted as 
the path the system follows in mesoscopic phase space 
in order to sustain a given current fiuctuation q. This 
path may be in general time-dependent, and the result- 
ing general variational problem is remarkably hard. This 
problem becomes simpler however in different limiting 
cases. For instance, one expects that small current fluc- 
tuations around the average, g ~ (g), result from the 
random superposition of weakly-correlated (if any) local 
fluctuations of the microscopic jump process. In this case 
one expects the optimal density fleld to be just the flat, 
steady-state one, pq{x,t) = po, and hence jq(x,t) — q, 
resulting in a simple quadratic form for the current LDF 



Gflat(9) 



{q-<T{po)EY 

Mpq) 



(3) 
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Therefore gaussian statistics is obtained for small (i.e. 
typical) current fluctuations, in agreement with the cen- 
tral limit theorem. The argument above breaks down 
however for moderate current fluctuations. In fact, Bod- 
ineau and Derrida have shown recently [TO that the flat 
profile indeed becomes unstable, in the sense that G{q) 
increases by adding a small time-dependent periodic per- 
turbation to the otherwise constant profile, whenever 
87T^D^{po)a{po) + (E^a^ipo) - q>"{po) < 0, where a" 
denotes the second derivative. This condition yields a 
critical current 



m = 



+ E^a^p,)., 



(4) 



which signals the onset of the instability. This instability 
can be interpreted as a dynamical phase transition at the 
fluctuating level, and involves the spontaneous breaking 
of translation symmetry (see Fig. [2]). In fact, for the 
WASEP the dynamic phase transition corresponds to the 
emergence of a macroscopic jammed state which hinders 
transport of particles to facilitate a current fluctuation 
well below the average. When the instability kicks in, an 
analysis of the resulting perturbation [1^ suggests that 
the dominant form of the optimal profile is a traveling 
wave, pq(x,t) = LOq{x — vt), moving at constant velocity 
V across the system [TBI HZ] ■ Provided that the traveling- 
wave form remains as the optimal solution for currents 
well-below the critical threshold, the current LDP can 
now be written as 



G{q) = 



dx 



Lj,ix),vJ„ 2a[u;q{x)] 



[q - vpo +VUJq{x) 



+D[u;q{xM{x) - a[uq{x)]E]\ (5) 



where the minimum is now taken over the traveling wave 
profile uig{x) and its velocity v. Notice that, for the insta- 
bility to exist, the strength of the driving field, \E\, must 
be large enough to guarantee a positive discriminant in 
eq. Q, namely 



\E\ > \Er, 



Re 



o-(Po)cr"(po) 



(6) 



Therefore, since the mobility <y{p) is positive definite, a 
non-zero threshold field only exists for models such that 
(t"(p) < 0, which is the case of the WASEP here stud- 
ied, where a{p) ~ p{l — p). Other transport models, as 
for instance the Kipnis-Marchioro-Presutti (KMP) model 
of heat conduction [17, 21j, have <j"{p) > and hence 
\Ec\ — 0, thus exhibiting the aforementioned instability 
even in the absence of external fields [17]. 

It is worth noting that MFT inherits the microscopic 
symmetries of the system of interest. In our particular 
case, the WASEP shows a clear particle-hole symmetry, 
and this is reflected in the current LDF, above and below 
the instability. In particular, the optimal wave profile 
ujq{x), associated with a current fluctuation \q\ < \qc\ 
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FIG. 2; (Color online) Typical evolution of microscopic con- 
figurations for current fluctuations above and below the crit- 
ical current for three different densities in the WASEP. Left 
panels correspond to currents above the critical one where 
the system remains homogeneous. Right panels correspond 
to subcritical current fluctuations where a traveling wave 
emerges. The velocity of the traveling wave of the top right 
panel (po = 0.3) is positive. The traveling wave of the cen- 
tral right panel does not move on average, corresponding to 
Po = 1/2, and the wave at the bottom left panel (po = 0.7) 
moves with negative velocity. 



for a density po, is complementary to the optimal wave 
profile for the same value of q and density 1 — po: i-e. 



Wq(x; Po) = 1 - i^qix; 1 - Po). 



(7) 



In addition, the optimal wave for density 1 — po trav- 
els with the same speed but opposite direction to that 
of the corresponding wave for density po, i.e. Wq(po) — 
— tig(l — Po). In the particular case of po — 1/2 these rela- 
tions imply that, for any current fluctuation, the optimal 
density profile and its complementary are equivalent, and 
the velocity of the optimal traveling wave is zero for all 
fluctuations. Therefore, for po = 1/2 the typical macro- 
scopic configurations in the time-dependent regime have 
a well defined wave structure which however does not 
move on average. This can be observed in Fig. [2| where 
typical system space-time trajectories for current fluctu- 
ations above and below the critical current are displayed 
for Po = 0.3, 0.5, and 0.7. Notice that for \q\ < \qc\ there 
is a nontrivial structure which travels with opposite ve- 
locities for Po = 0.3 and 0.7, and which does not move 
when Po = 0.5. Furthermore, using the above symmetry 
relations for the WASEP current LDF we find that 



G(g;po) = G{q; 1 - po) ■ 



(8) 



Hence, given a external field, it is enough to compute the 
current LDF for po G [0, 1/2]. 

Another interesting symmetry, though far less obvious, 
is related to the time-reversibility of microscopic dynam- 
ics. This relation, known as the Gallavotti-Cohen flue- 
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FIG. 3: (Color online) Legendre transform of the current 
LDF, /i(A). Top: Measured fJ,{X) for po = 0.3 and increas- 
ing TV, together with the MET result (solid red line) and the 
quadratic approximation (dashed blue line). Bottom: Equiv- 
alent data for po = 1/2. Insets: p(A) — /xflat(A) for the same 
A*' and po — 0.3 (top) and 0.5 (bottom). In all cases, data 
converge to the MET prediction as A'" increases. 



tuation theorem [P-^, implies a remarkably simple con- 
nection between the probability of a given current fluc- 
tuation q and the reverse event, —q, which can be stated 
for the current LDF in the following way 



III. NUMERICAL RESULTS 

Our aim in this paper is to characterize in detail the 
dynamical phase transition in current statistics for the 
WASEP using numerical simulations, and compare with 
the predictions derived within MFT. These predictions 
are explicitly worked out in the Appendix. The critical 
current and the threshold field for the WASEP are 



VPo(l - Po) 



(11) 

(12) 



Equivalently, taking into account that (q) — po{l — po)E, 
we have that 



l(9>K 1 



(13) 



Giq)-Gi-q) = 2Eq. 



(9) 



SO \qc\ < 1(g) I in WASEP for field strengths above the 
threshold, while no phase transition happens for < 
\Ec\. 

In order to investigate the instability described in Sec- 
tion [IT] using numerical simulations, we need to explore 
the statistics of both typical and rare current fluctua- 
tions. While the former pose no problem and can be 
studied in standard simulations, to sample the atypical 
trajectories associated with rare current fluctuations we 
must resort to advanced Monte Carlo methods that allow 
to measure directly LDFs in many particle systems [221- 
[24] . This technique implies a modification of the stochas- 
tic microscopic dynamics, in such a way that the rare 
events responsible of a large current fluctuation are no 
longer rare with the modified dynamics. The numerical 
method also requires the parallel simulation of multiple 
clones or copies of the system which may be 

replicated or pruned depending on its importance for the 
particular fluctuation we want to measure. In this work 
we used in particular Nc — 2 x 10^ clones for pq = 0.3 
and Nc ^ 5 X 10^ for po = 1/2, and we checked that 
results do not depend on the total number of clones for 
large enough Nc [25]. The method yields a Monte Carlo 
estimate of the Legendre transform of the current LDF 



This in turn implies that the odd part of the typically 
nontrivial function G{q) is linear in the current, with 
an universal coefficient 2E. This symmetry can be also 
stated for the Legendre transform of the current LDF, 
see eq. ( 14 ) below 



p{X) =p{-X-2E). 



(10) 



As we will see below, this fluctuation relation is fully 
confirmed in our simulations, both below and above the 
instability. Moreover, the Gallavotti-Cohen relation can 
be used to bound the validity of the simulation method 
used to explore large deviations [131 [5S] , see below. 



p{X) = TaBx[G{q) + \q\ , 

9 



(14) 



with A a parameter conjugated to the current, such that 
G'{q) -I- A = 0. In this way, the function ^(A) can be seen 
as the conjugate potential to G{q), a relation equivalent 
to the free energy being the Legendre transform of the 
internal energy in thermodynamics, with the tempera- 
ture as conjugate parameter to the entropy. Legendre- 
transforming the quadratic current LDF in eq. ([s]) 
once particularized for WASEP, obtained for the time- 
independent (homogeneous) fluctuation regime, we have 

Mflat(A) = ^(l-po)A(A + 2£;). (15) 
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threshold Ec for all three densities po, see eq. (12), so 



FIG. 4: (Color online) Top: Large Deviation Function for 
Po = 0.3 Inset: Measured average current qx as a function of 
A and increasing N, together with the analytical prediction 
base on the MFT. Bottom: Same results for po = 1/2. For 
both densities the traveling wave solution enhances the proba- 
bility for fluctuations \q\ < \qc\ (solid red line) with respect to 
the flat profile associated to gaussian statistics (dashed black 
line) . 



< Wc 



However, for currents in a well-defined interval \q\ 
with \qc\ defined in eqs. (Ill or (13), a time-dependent 



regime with an optimal density profile in the form of a 
traveling wave is expected. This corresponds to values of 
the conjugate parameter A such that |A -I- i?| < Ac, with 



Po(l - Po) 



Po(l - Pq) 



(16) 

and where we have used eqs. (11|-(13). Therefore we 



expect traveling wave solutions for < A < A+, where 
Xf=±A,,-E. (17) 
In this way the time-dependent fluctuation r egim e kicks 



in whenever /Li(A) < p-{X^) — —'k /2, see eq. (151 



We performed simulations of the Id periodic WASEP 
for three different average densities, po — 0.3, 0.5 and 
0.7, for increasing system sizes N G [8,64] and a fixed 
external field E = +10. This driving field is above the 



we expect the instability to appear on the basis of the 
analysis of Section |IT] Fig. [3] shows simulation results 
for /i(A) and increasing values of N for two different val- 
ues of Pq, together with the explicit MFT results (see 
the Appendix). Gaussian current statistics corresponds 
to a quadratic behavior in p{X) paatW, see eq. (151, 
which is fully confirmed in Fig. [3] for |A -I- i?| > Ac and 
different values of N. This means that small current fluc- 
tuations, g ~ (9), have their origin in the superposition of 
uncorrelated (or at most weakly-correlated) local events 
of the stochastic jump process, giving rise to Gaussian 
statistics as dictated by the central limit theorem and 
thus confirming their incoherent origin. Interestingly, 
this observation also applies to the time-reversal part- 
ners of these small fluctuations, —q, which are far from 
typical. This is a direct consequence of the Gallavotti- 
Cohen symmetry [IHl], which implies that the statistics 
associated with a current fluctuation does not depend on 
the current sign [5S]. On the other hand, for fluctua- 
tions below a critical threshold, |g| < \qc\ or equivalently 
A^ < A < Ajf , deviations from this simple quadratic form 
are apparent, signaling the onset of the dynamical phase 
transition anticipated in Section [llj see also Fig. [2] In 
fact, as N increases a clear convergence toward the MFT 
prediction (which is strongly non-quadratic in the regime 
A~ < A < A+) is observed, with very good results already 
for N — 64. Strong finite size effects associated with the 
finite population of clones Nc prevent us from reaching 
larger system sizes [131 HZl US] , but iV = 64 is already 
close enough to the asymptotic hydrodynamic behavior. 
Still, small corrections to the MFT predictions are ob- 
served, see the insets of Fig. [3j which quickly decrease 
with N. On the other hand, the Gallavotti-Cohen fluc- 
tuation theorem for currents holds in the whole current 
range, see eq. (10 1 and Fig. [sj both in the homoge- 
neous and time-dependent current fluctuation regimes. 
Furthermore, the Gallavotti-Cohen symmetry holds irre- 
spective of N, as a result of the microreversibility of the 
model at hand: while we need a large size limit in order 
to verify the predictions derived from MFT (which is a 
macroscopic theory), no finite-size corrections affect the 
fluctuation theorem, whose validity can be used to ascer- 
tain the range of applicability of the cloning algorithm 
used to sample large deviations pS] . 

We also measured the average current (qx) associated 
with a fixed value of the conjugate parameter A. The in- 
sets of Fig. |4] show our results and the predictions based 
on MFT. Again the agreement is excellent (improving as 
N increases) , both above and below the dynamical phase 
transition, even though there is a clear change of behav- 
ior across the transition points Xf. In particular, in the 
Gaussian fluctuation regime the current is linear in A, 
namely {q\) ~ po(l ^ Po)i^ + E), while the relation be- 
comes strongly non-linear in the time-dependent region, 
|A + i?| < Ac. We may use the measured (qx) to give 
a direct Monte Carlo estimate of the current LDF G{q). 
In fact, (qx) is the current conjugated to a given A and 




FIG. 5: (Color online) Measured traveling wave profiles for different values of A G (AJ , A J), varying A'^ £ [8, 64] and two different 
average densities po = 0.3, 0.5, together with theoretical predictions from MFT, see the Appendix, (a) uj\{x) measured for 
po ~ 0.3, different A and increasing A'^, and MFT predictions, (b) Measured density profiles as a function of A for A*' — 64. 
Optimal profiles are flat up to a critical current (equivalently A^) where a traveling wave emerges. The inset shows the MFT 
prediction, (c) and (d) are equivalent to (a) and (b), respectively, but for an average density po = 0.5. Notice the comparatively 
thicker wave for po = 0.5. (e) Collapse of measured profiles associated with different current fluctuations (jJ\{x) and their time- 
reversal partners uJ-x-2e{x) for A'^ = 64 and po = 0.3, together with theoretical predictions. Optimal profiles, both below and 
above the instability, remain invariant under change of sign of the current, (f) Same results for po — 0.5. 



hence we may write G{q) — /i(A) + X{q\), where we com- 
bine the measured /i(A) in Fig. [sjand the measured (^a) 
in the insets of Fig. |4 The result for G{q) and different 
values of po is plotted in Fig. |4j where we again find a 
good agreement between theory and Monte Carlo simu- 
lation. Notice in particular the deviation from quadratic 
behavior observed for currents \q\ < \qc\ resulting from 
the formation of macroscopic jammed states (see below). 

The dynamical phase transition is most evident at the 
configurational level, as observed in Fig. [2] so we mea- 
sured the average density profile associated with a given 
current fluctuation pJ5^ for different values of the total 
density po] see Fig. [5] Because of the system periodic- 
ity, and in order not to blur away the possible structure 
present in microscopic configurations, we performed pro- 
file averages around the instantaneous center of mass. 
For that, we consider the system as a Id ring embedded 
in two-dimensional space, see Fig. [l]a, and compute the 
angular position of the center of mass, shifting it to the 
origin before averaging. In particular, we assign an an- 



gular position 9i = 2m /N to each site i € [1, A'^] in the 
lattice. The angular position of the center of mass for a 
given microscopic configuration n = {ui^i = 1, . . . , A'^}, 
with Ui = 0, 1 the on-site occupation numbers, is thus 
defined as 



M 



tan- ('^) 



(18) 



where 



1 ^ 

XcM = — cos ( 

1=1 



CM 



1 ^ 

M ^ 

i=l 



Tij sin 9j 



(19) 

and recall that M = X]i=i the total number of parti- 
cles. Notice that this center-of-mass averaging procedure 
yields a spurious weak structure in the Gaussian (ho- 
mogeneous) fiuctuation region, equivalent to averaging 
random particle profiles around their (random) center of 
mass. Such a spurious profile is of course independent 
of the current q and can be easily subtracted. On the 
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FIG. 6: (Color online) Measured velocity as a function of A 
for po = 0.3, po = 1/2 and increasing N, together with the 
MFT result. 



Fig. [2] and making statistics for the measured velocity. 
Figure [6] shows the mean velocity for At = 200 Monte 
Carlo steps as a function of A for increasing values of 
iV and different values of poi ^nd the agreement with 
MFT is again very good already for = 64 (we checked 
that other values of At yield equally good results). As 
for the current, see the insets of Fig. |4j there is a clear 
change of tendency across the instability points A J , with 
the velocity being a linear function of A in the Gaussian 
regime but turning strongly nonlinear in the traveling- 
wave region. Interestingly, the measured wave average 
velocity for po = 0-5 is compatible zero for all current 
fluctuations, hence confirming the MFT prediction based 
on the particle-hole symmetry of WASEP. 



IV. CONCLUSIONS 



other hand, once the instability is triggered average pro- 
files exhibit a much more pronounced structure resulting 
from the appearance of a traveling wave; see right col- 
umn in Fig. [2j Figures [5ja,c show the measured profiles 
ij:\{x) for different A G (A7,A+), varying N G [8,64] 
and two different average densities po = 0.3 (a) and 0.5 
(c), together with theoretical predictions from MFT as 
calculated in the Appendix. Again, fast convergence to- 
ward the MFT result is observed, with excellent agree- 
ment for = 64 in all cases. Moreover, Figs. [5]b,d 
show a three-dimensional plot of the measured profiles 
for TV = 64 and different A, again for po = 0.3 (b) and 
0.5 (d), which closely resembles the MFT scenario plotted 
as insets to these figures. In general, the traveling wave 
profile grows from the flat form as A crosses the critical 
values A^ penetrating into the critical region, thus favor- 
ing a macroscopic jammed state that hinders transport 
of particles and thus facilitates a time-averaged current 
fluctuation well below the average current. This macro- 
scopic jammed state reaches its maximum expression for 
A = —E, or equivalently (7 = -see the insets of Fig. 
[4j so the system is maximally jammed for zero current 
irrespective of the driving field E. 

An interesting corollary of the Gallavotti-Cohen fiuc- 
tuation symmetry [THl] is that the optimal density 
profile associated with a given current fluctuation re- 
mains invariant under changes of the current sign, i.e. 
Lt^qix) = LO-q{x)^ independently of the driving external 
field [ini US]. We confirm this property in Fig. [5]e,f, 
where we plot for = 64 and different po the optimal 
density profiles for different pairs of fluctuations coupled 
by time-reversal, i.e. for pairs of values (A, — A— 2i?), find- 
ing an excellent collapse as predicted by theory. More- 
over, the collapsing pairs of profiles agree to a high degree 
of accuracy with the theoretical curves. 

We also measured the average velocity associated with 
a given current fluctuation by fltting the motion of the 
center of mass during small time intervals At to a ballistic 
law, rcin(t -I- At) — rcm(i) = vt^ see, e.g., right column in 



In this paper we have considered the statistics of the 
time-averaged current in a simple driven diffusive system, 
the Id weakly asymmetric simple exclusion process on a 
ring. This system exhibits a dynamical phase transition 
at the fluctuation level for large enough driving external 
fields, and we have characterized such instability in detail 
using advanced numerical simulations. We show in par- 
ticular that typical (i.e. small) current fluctuations result 
from the sum of weakly-correlated local events which take 
place in an otherwise homogeneous density background 
on average, and thus give rise to Gaussian current statis- 
tics in agreement with the central limit theorem. How- 
ever, for large enough current fluctuations well below the 
average current (namely in a well-deflned range of current 
fluctuations around zero current), the system breaks this 
homogeneity and self-organizes its density proflle into 
a macroscopic jammed state in the form of a coherent 
traveling wave moving at constant velocity. This wave 
structure hinders transport of particles and thus facili- 
tates a time-averaged current fluctuation well below the 
average current. It is worth emphasizing that the emer- 
gence of a traveling wave breaks spontaneously a sym- 
metry (translation invariance) at the fluctuating level. 
This phenomenon is fully captured by macroscopic fluc- 
tuation theory, whose predictions are explicitly worked 
out and confirmed to a high degree of accuracy by sim- 
ulation results. Our study offers insights not only on 
the current large deviation function which characterizes 
current statistics, but also on the optimal density profiles 
associated with the different fluctuations and their veloc- 
ity, as well as the effect of flnite-size corrections on the 
observables of interest, providing a detailed characteriza- 
tion of the dynamical transition not available before. 

A similar dynamical phase transition has been re- 
cently observed and characterized in another model of 
transport, the Kipnis-Marchioro-Pressuti (KMP) model 
of heat conduction [17] . In that case the phenomenon is 
even more striking, as it happens in an isolated system in 
the absence of any external fleld, spontaneously breaking 
a symmetry in ID and illustrating the idea that criti- 
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cal phenomena not allowed in equilibrium steady states 
may however arise in their fluctuating behavior or under 
nonequilibrium conditions. Although both instabilities 
are described equally well by MFT, the physical inter- 
pretation of the dynamical transition is quite different. 
In particular, in the KMP model the instability hap- 
pens because the system optimizes the transport of a 
large current by gathering energy in a localized packet 
(the wave) which then travels coherently, breaking spon- 
taneously translation symmetry in the process. On the 
other hand, the WASEP instability happens in order to 
hinder the transport of particles via the formation of a 
macroscopic jammed states (the wave), thus facilitating a 
current fluctuation well below the average. Interestingly, 
both phenomena are sides of essentially the same insta- 
bility. It is also worth noticing that similar instabilities 
have been described in quantum systems j27j . 



Our results show unambiguously that the dynamical 
phase transition observed in the WASEP is continuous 
as conjectured in [T^ , excluding the possibility of a first- 
order scenario, in concordance with previous observations 
for the KMP model ^ITj. This suggests that a travel- 
ing wave is in fact the most favorable time-dependent 
profile once the instability is triggered. This observa- 
tion may greatly simplify general time-dependent calcu- 
lations, but the question remains of whether this is the 
whole story or if other, more complex solutions may play 
a dominant role for even larger fluctuations. An inter- 
esting, related question concerns the properties of time- 
dependent solutions for systems with open boundaries, 
where traveling-wave patterns are not appropriate. The 
time-independent profiles in these cases, from which a 
suitable perturbation analysis would hint at the form of 
the time-dependent solution, are far more complex than 
the trivial homogeneous profiles that appear for periodic 
systems, difficulting progress along this line. In fact, a re- 
cent study p8i has found no evidence of dynamical phase 
transition in WASEP with open boundaries. In any case, 
it seems clear that extremely rare events call in general 
for coherent, self-organized patterns in order to be sus- 
tained l29l. 



Another interesting direction to explore in a near fu- 
ture is the appearance of this phenomenon in higher- 
dimensional systems. In this case the solution of the 
associated MFT is far more complicated, with no guar- 
antee of an unique solution and several competing pat- 
terns already known |30) . The role of numerical simu- 
lation will hence prove essential to explore rare current 
fluctuations in high-dimensional systems and to under- 
stand the appearance of dynamical phase transitions at 
the fluctuation level [30] • Furthermore, the simplicity 
and elegance of this phenomenon suggests that it might 
be a rather general property of any fluctuating fleld the- 
ory, with possible expressions in quantum field theory, 
hydrodynamics, etc. 
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Appendix: Macroscopic fluctuation theory for 
current statistics and dynamical phase transition 

Macroscopic fiuctuation theory (MFT) |10 describes 
in detail dynamic fluctuations in driven diffusive sys- 
tems, starting from the hydrodynamic evolution equa- 
tion for the system of interest and the sole knowledge 
of two transport coefficients, which can be measured ex- 
perimentally. From this knowledge, MFT offers explicit 
predictions for the current LDF and the associated path 
in phase space responsible of a given fluctuation. MFT 
applies to systems described at the mesoscopic level by a 
(fluctuating) continuity equation of the form 



dtp + dxj = , 



(A.l) 



where p{x, t) and j(a;, t) are the density and current flelds, 
respectively, and t and x £ [0, 1] are the macroscopic 
time and space variables, obtained after a diffusive scal- 
ing limit such that x — i/N and t = i/N'^, with i and 
t the microscopic space and time variables. Periodic 
boundary conditions, the case of interest here, thus im- 
ply p{0, t) = p{l, t) and j(0, t) = t). Moreover, as the 
system is isolated, the total density remains constant 



Po 



p{x, t)dx . 



(A.2) 



The current fleld in eq. (A.l I is in general a fluctuating 
quantity, and can be written as 

j{x, t) = ~D{p)d^,p{x, t) + a{p)E + e(x, t). (A.3) 

The flrst term is Pick's (or equivalently Fourier's) law, 
where D{p) is the diffusivity (which might be a nonlinear 
function of the local density). The second term is just 
the coupling to the external fleld E, mediated by the so- 
called mobility <j{p), and S,{x,t) is the current noise that 
is gaussian and white. 



{ax,t)) = o, {ax,mix\t')) 



N 



S{x-x')6{t-t'), 
(A.4) 

This gaussian fluctuating field is expected to emerge for 
most situations in the appropriate mesoscopic limit as a 
result of a central limit theorem: although microscopic 
interactions for a given model can be highly complicated, 
the ensuing fluctuations of the slow hydrodynamic fields 
result from the sum of an enormous amount of ran- 
dom events at the microscale which give rise to gaus- 
sian statistics, with an amplitude of the order of N~-^/^, 
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in the mesoscopic regime in which eq. (A.l) emerges 



For long times, a system described by the above set of 
equations reaches a nonequihbrium steady state charac- 
terized by a homogeneous density distribution po a-nd a 
nonzero net average current (q) — a[po)E. Note that, 
for the WASEP, the two essential transport coefficients 
D{p) and cr(p), which determine the complete macro- 
scopic fluctuating behavior of the system, are D{p) = ^ 
and a{p) = p{l — p) [H]. In what follows we describe the 
theory in general, only particularizing for the WASEP 
case in the last stages of the calculation. 

A simple path integral calculation starting from eq. 
(A.l) then shows that the probability of a given history 



or path in mesoscopic phase space (i.e. the space 

spanned by the hydrodynamic fields) obeys a large devia- 
tion principle of the form j}o) ~ exp(— A^I^[p, j]), 
where the rate function is given by [lOl [131 HI! 



Mil) 



Jo Jq 



We are interested here in the fiuctuations of the space 
and time integrated current 



dt / dxj{x, t) . 

^0 



(A.6) 



The probability of observing a given q can now be writ- 
ten as a path integral over all possible histories {p,j}o, 
weighted by its probability measure Pi{j, p}q), and re- 
stricted to tho se his torie s com patible with the value of q 
and po in eqs. (A.6) and ( A.2), respectively, and the con- 



tinuity equation (A.l) at every point of space and time. 



For long times and large system sizes, this path integral 
is dominated by the associated saddle point and scales as 
P{q) ^ exp{+TNG{q)}, where G{q) is the current large 
deviation function (LDF) given by 



G(g) = - lim 



- min (p, j) 

T {pj}S 



(A.7) 



The optimal density and current fields solution of this 
variational problem, pq(x,t) and jq{x,t), can be inter- 
preted as the optimal path the system follows in order to 
sustain a long-time current fluctuation. Finding the opti- 
mal fields is in general a complex spatiotemporal problem 
whose solution remains challenging in most cases. The 
problem becomes much simpler however in different lim- 
iting cases. For instance, one expects that small current 
fluctuations around the average, q ~ {q), result from the 
random superposition of weakly-correlated local fluctu- 
ations of the microscopic jump process. In this case it 
is reasonable to assume the optimal density fleld to be 
just the flat, steady-state one, Pq{x,t) = pq, and hence 
jq{x,t) — q, resulting in a simple quadratic form for the 
current LDF 



Gfiat(g) 



{q'<Tipo)Ey 

Mpq) 



(A., 



Therefore Gaussian statistics is obtained for small (i.e. 
typical) current fluctuations, in agreement with the cen- 
tral limit theorem. The previous argument, revolving 
around small fluctuations, breaks down however for mod- 
erate current deviations where correlations may play a 
relevant role. In fact, Bodineau and Derrida have shown 
recently [T5] that the flat proflle indeed becomes unsta- 
ble, in the sense that G{q) increases by adding a small 
time-dependent periodic perturbation to the otherwise 
constant profile, whenever 

8Ti''D\pa)a{po) + (E^a^ipo) - q^ViPo) < , (A.9) 

where cr" denotes second derivative. This condition im- 
plies a well-defined critical current 



I S7T^D^Po)cTipo) 
^"{Po) 



E^a^po)^ 



(A.IO) 



In the WASEP case this imphes that for any \q\ < jgd the 
instability emerges. This instability can be interpreted as 
a dynamical phase transition at the fluctuation level, and 
involves the spontaneous breaking of translation symme- 
try (see Fig. [2]). In fact, the formation of a traveling wave 
corresponds to the emergence of a macroscopic jammed 
state which hinders transport of particles to facilitate a 
current fluctuation well below the average. Notice that, 
for the instability to exist, the strength of the driving 
fleld, \E\, must be large enough to guarantee a positive 
discriminant in eq. (A.IO), namely 



\E\ > 1^,1 = Re 



a{po)(T"{po) 



(A.ll) 



Therefore, since the mobility <y{p) is positive deflnite, a 
non-zero threshold field only exists for models such that 
ct"(p) < 0, which is the case of the WASEP here stud- 
ied, where cr(p) = p{l — p). Other transport models, as 
for instance the Kipnis-Marchioro-Presutti (KMP) model 
of heat conduction [TTlini], have cr"{p) > and hence 
\Ec\ — 0, thus exhibiting the aforementioned instability 
even in the absence of external fields [XT' . 

When the instability kicks in, an analysis of the result- 
ing perturbation [TQ suggests that the dominant form of 
the optimal profile is a traveling wave moving at constant 
velocity v 



Pq{x,t) = UJq{x - Vt) , 



which implies via the continuity equation (A.l) 

jq{x, t) = q - Vpo + VUJq{x - vt) . 



(A.12) 



(A.13) 



Provided that the traveling- wave form remains as the op- 
timal solution for currents well-below the critical thresh- 
old, the current LDF can now be written as 



G{q) 



dx 



[q - vpo + VUJq{x) 



+D[ujq{x)]u:'q{x) - a[LOq{x)]Ef (A.14) 
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where we have dropped the time dependence due to the 
periodic boundary conditions, and the minimum is now 
taken over the traveling wave profile uJi, {x) an d its veloc- 
ity V. Expanding now the square in eq. (A. 14), we notice 
that the terms linear in give a null contribution due 
again to the system periodicity. Taking also i nto account 



the constraint (jjq{x)dx — po, see eq. (A. 2), one gets 



Giq) 



mm 



dx{X{ujg)+uj'{xfY{ujq)) 



where, borrowing the notation of ref . [16| , 

[q ~ v{po - ujg)]^ , E^a{ujg) 



and 



2a{ujq) 



(A.16) 



(A.17) 



The differential equation for the optimal profile solution 



of the variational problem eq. (A. 15 1 can be written as 

Xiijq) ~ u;'g{xfY{u;g) = Ci + . (A.18) 

This equation gencrically yields a symmetric optimal pro- 
file with a ujq(x) with a single minimum uji — ijjq{xi) and 
a single maximum luq — u!q{xQ) such that |a;o — xi \ = 1/2 
pS] . The constants Ci and C2 can be expressed in terms 
of the extrema wi and wq. 

The optimal velocity also follows from the above vari- 
ational problem, 



dx 



(j{ujq) 



dx 



{Ujq - PqY 
(j{ujq) 



(A.19) 



It is worth emphasizing that the optimal velocity is pro- 
portional to q. This implies that the optimal profile 
solution of eq. (A.18 1 depends exclusively on and 



not on the current sign, reflecting the Gallavotti-Cohen 
time-reversal symmetry. This invariance of the optimal 
profile under the transformation q -H- —q can now be 



used in eq. (A. 15) to show explicitly the GC symmetry 
G{q) — G{—q) = 2Eq. This fluctuation relation is fully 
confirmed in the simulations discussed in the main text. 



The constants Ci and C2 appearing in eq. (A.18) can 



be expressed in terms of the extrema Wi and of the 
profile via 



X{ijji) = Ci+ , 

a: (wo) = Ci + C2L00 . 



(A.20) 
(A.21) 



Moreover, the extrema locations are fixed by the con- 
straints on the distance between them and the total den- 
sity of the system. 



dx 




Y{^q) 



X{ujq) - Cl - C2UJq 



dujq (A.22) 



and 

Po 
2 



UJq{x)dx 




^lY{u:q) 



-duj„ 



C2Wa " 



X{Uq) - Cl 

(A.23) 

where we have used in the last equality of both expres- 
sions the differential equation (A.18). In this way, for 



fixed values of the current q and t he den sity po (provided 
qE, externally), we use eqs. ( A.19)-(A.23) in order to de- 
termine the five constants wi, cjq, Ci, C2, w which can be 
(A. 15) used in turn to obtain the shape of the optimal density 



profile ujq{x) from eq. (A.18). 



Notice that the unknown variables ujo, uji appear as 
the integration limits in eqs. (A.22) and (A.23), making 



this problem remarkably difficult to solve numerically. 
In what follows we show how, by performing a suitable 
change of variables, the integrals involved in the calcu- 
lation can be transformed into known functions, as e.g. 
elliptic integrals of the first kind, thus allowing to derive 
an explicit analytical expression for ujq{x) as a function 
of the relevant constants. We start by doing a change 
of variables to express all the relevant magnitudes in di- 
mensionless form 



Po 



u; ujq{x) = poh{x) : 



E: 



1 " 

Po 



(A.24) 



Particularizing now our calculation for the WASEP, 
where the transport coefficients are D(p) = 1/2 and 



cr(p) = p(l — p), the differential equation (A.18) for the 
traveling wave reads 

h!{x) = ^ [(1 - u + uhf - 2D2h^{l - poh)~ 
Po 



2Dih{l - poh) + -^h^{l - pohf 
Po 



1/2 



where we have defined -D, 



(A.25) 
1, 2. Moreover, 



eqs. ( A.20 )-( A.21 ) can now be written as 
(1- 



■ uhk) 



2D2hi{l - pohk) + 2D^hk{l - pohk) 

^hlil-pohk)\ (A.26) 
Po 



where h^. = loi^/pq, with fc = 0, 1, are the extrema of the 
dimcnsionless profile h{x). We can now use the above 
equations to write the constants Di and D2 as a function 
of the dimcnsionless variables hi, ho, e and u, 



D2 = 



1 - u 

2(1 - poho){l - pohi 
(1-u) 



(1 -u 
hiho 



[1 - Po{ho 



hi)] 
(A.27) 



e 



.[1 



Poiho + hi)], 



[1 — u -\- uho)"^ 
2/io(l - poft-o) 



-D2h 



2"-0- 



feo(l - pofep) 
pI 2 



(A.28) 
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The remaining task consists in obtaining the three un- 
known variables hi, ho and u from eqs. (A. 19 1, (A. 22) 



and (A.23). In particular, eq. (A.22) boils down to 



dh[{l - u + uh)^ - 21)2/1^(1 - poh)- 



hi 



2DM^ - Poh) + -^h\l - pohf]-'/^ 
Po 



Po' 



(A.29) 



The integrand of the above expression can be written in 
the following product form 



{-ah + b){h + c)ih~ hi){ho - h), 



(A.30) 



where the coefficients a, b and c (obtained by matching 
order by order) are (simple) functions of the unknown hi, 
ho, u and the known q, po and e. In order to eliminate the 
unknown extrema hi, ho from the integration limits in 
eq. (A.29), we perform the following change of variables 



h = ho — a{ho — hi). 



(A.31) 



which allows us to rewrite eq. (A.29) as 

»i 



Po 



da 

{ho~hi)y/a 



h 







ho - hi 



X (1 — a)a a 



aho — b 
a{ho - hi) 



-1/2 



Defining now 

2 _ (ac + b){ho - hi) 



{~ahi+b){ho+c) 



and z 



ho - hi 
772/11 



(A.32) 



(A.33) 



we get that Eq. (A.32) turns into 

q _ 2 da 

Po y/{ac + b)zhi Jo y^{l - a^){l - rfa'^) 
2 



{ac + b)zh 



(A.34) 



where K{rf) is the complete elliptic integral of the first 
kind. It is worth emphasizing that a, b, c, z, and 77^ 
depend on hi, ho, u and on q, po and e. 

In a similar way, we can derive an expression for the 
adimensional optimal profile h{x) by writing 



dx — X ~ [ ^z— 
Jhi h' 



(A.35) 



and proceeding in the same way as before. This yields 



qx 



1 



da 



Po ^/{ac + b)zhl J J y(I^^Q;2y(j"Tr^^2;^ 

which is equivalent to 

qx K{ri'^) — F[sin^"^(7), 77^] 
Po 



(A.36) 



where 



■y/ {ac + b)zhi 



l{x) 



^ {ho + c)[h{x)-hi] 
{ho - hi)[h{x)+c] 



(A.38) 



and 



^[sin-i(7),772] = 



da 



V(l 



(A.39) 



is the incomplete elliptic integral of the first kind. Now, 
by using eqs. (A.34) and (A.37) we deduce that 



F[s\n-\-i),r^^]=K{r^^){l-2x). 



(A.40) 



Solving for 7 we obtain 

7(0;) = JacobiSN [i^(?7^)(l - 2x)] , (A.41) 



where JacobiSN is the inverse of the incomplete elliptic 
integral of the first kind. This yields finally the optimal 
density profile, uj{x) = poh{x), with h{x) obtained from 
the above equation after taking into account eq. (A.38) 



h{x) 



hi + cT{x) 
1 - T(a;) 



(A.42) 



with T{x) = (JacobiSN [i^(r7^)(l - 22;)]) 



2 ho — hi 
ho + c 



Equation (A.42 ) for the optimal traveling wave reflects 



the explicit dependence of the wave proflle on the con- 
stants hi, ho and u. The remaining job consists in solv- 
ing numerically for these constants in a self-consistent 
manner, once the explicit dependence of the extrema has 
been removed from integral limits. These constants can 
be thus obtained from eqs. ( |A.19[ ), ( |A.22[ ) and ( |A.23[ ), 
which can be written as 



dx^^^^^K^{l ~u + uh{x)) - 0, 



h{x)^ 



^ {ac + b)hiz 



K{v') 



<]_ 
Po'' 



h{x)dx — 



2' 



(A.43) 



(A.44) 



(A.45) 



where a, b, c, 77^ and z are known functions of hi, ho, u. 
In this way, for given values of q, po and E, we get self- 
consistently hi, ho, u using the form of the profile ob- 



tained in eq. (A.42), which depends explicitly on these 
constants. 

To obtain the current LDF, G{q), we just integrate 



numerically its expression (A.15) once particularized for 
the WASEP (cr(a;) = a;(l - w)7~D(a;) = \), using the op- 
timal wave profile and velocity obtained from the previ- 
ous calculation. Finally, to compute the Legendre trans- 
form of the current LDF, we just evaluate numerically 
^(A) = max^[Ag -I- G{q)] = \q* + G{q*) with q*{X) solu- 



(A.37) tion of the following equation 



dG{q) 



dq 



q=q* 



\^ q*{X)-v{po-^{x)) ^ 
uj{x){l - uj{x j) 

(A.46) 
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